function [] = runParamRecovery_generate(runId, NumReplications)
rng(1); % for replicability

%%% Add paths
addpath('../general_funcs');
addpath('../model_funcs');
addpath('../model_funcs/dim_specific_funcs/dim2');
addpath('../model_funcs/dim_specific_funcs/dim3');
addpath('../consumerAdoptions');
addpath('../providerAdoptions');
addpath('../providerExits');
%%% Get input / output paths
config = getConfig_jointProcess(runId);

% Get config_consumerAdoptions and config_providerAdoptions: necessary input for make_Xdynamic_xxx functions
indivConfigs.consumerAdoptions = getConfig_consumerAdoptions(config.runId_consumerAdoptions);
indivConfigs.providerAdoptions  = getConfig_providerAdoptions(config.runId_providerAdoptions);
indivConfigs.providerExits      = getConfig_providerExits(config.runId_providerExits);

indivRunFolders.providerAdoptions  = config.runFolder_providerAdoptions;
indivRunFolders.providerExits      = config.runFolder_providerExits;
indivRunFolders.consumerAdoptions = config.runFolder_consumerAdoptions;

outputFolder  = createFolderIfNotExist(sprintf('%s/paramRecovery/continuous_time', config.runFolder)); % will store continuous-time data
outputFolder2 = createFolderIfNotExist(sprintf('%s/paramRecovery/day_level', config.runFolder)); % will store data aggregated at daily level

%%%%% Generate data: take as basis the initial conditions and observed covariates from the runDdata of the 3 models
% (of course the covariates that depend on installed bases will be determined endogenously)
%%% Load data from the 3 submodels
data0 = load_data_fullmodel(indivRunFolders);

%%% Get initial conditions (from base case) + Xdynamic_staticInputs
[initialConditions, Xdynamic_staticInputs0] = load_base_initial_conditions_and_Xdynamic_staticInputs(indivRunFolders);

%%%  Load parameter estimates from the 3 submodels => these are going to be the "true parameters" that generate Y in the synthetic data
[truth.params,paramsParts] = load_param_estimates(data0, indivRunFolders);

% If data0 used in runId is in continuous time, "aggregate it" at the day level to get non-network covariates at day level
if isfield(data0.consumerAdoptions, 'spellDurations')
	%%% Aggregate the data from spell level to day level
	data0.consumerAdoptions = aggregateCtsTimeData_consumerAdoptions(data0.consumerAdoptions);
	data0.providerAdoptions  = aggregateCtsTimeData_providerAdoptions(data0.providerAdoptions);
	data0.providerExits      = aggregateCtsTimeData_providerExits(data0.providerExits);
end

%%% Read dimensions
T  = data0.consumerAdoptions.dims.dims1{4};
K1 = data0.consumerAdoptions.dims.dims1{5};
K2 = data0.consumerAdoptions.dims.dims1{6};
K  = data0.providerAdoptions.dims.dims1{3};

surfaces = data0.surfaces;


%%% Combine dims into object mydims
mydims.T  = T;
mydims.K  = K;
mydims.K1 = K1;
mydims.K2 = K2;

%%% Compute the parts of V that does not depend on installed bases
XbetaFixed = calc_Xbeta_fixed(data0, paramsParts);

%%% Preprocess Xdynamic_staticInputs0 inputs to compute the parts of V that depend on installed bases
Xdynamic_staticInputs2 = permute_Xparts_static_X_X_FEs(Xdynamic_staticInputs0, mydims);


%%% Loop over replications
tic
for iter = 1:NumReplications
	disp(sprintf('Starting replication %d', iter));
	startingRandomState = rng;
	
	Xdynamic_staticInputs = Xdynamic_staticInputs0;
	
	%%%%% Generate Y based on continuous-time sequence and store timestamps
	[Y_full,~,spellDurations, spell2period] = generate_fullmodel_sequence_ctsTime2(indivConfigs, initialConditions, XbetaFixed, Xdynamic_staticInputs2, paramsParts);
	NumSpells = length(spellDurations);
	
	%%% Reconstruct installedBases and popAtRisk that were basically generated although not stored (using initialConditions and Yfull)
	disp('Reconstructing installedBases and popAtRisk...');
	% Make laggedProviders and popAtRisk for providerAdoptions and providerExits
	providerDeltas = Y_full.providerAdoptions - Y_full.providerExits; % NumSpells x K
	installedBases.laggedProviders = initialConditions.laggedProviders' + cumsum(providerDeltas, 1) - providerDeltas; % NumSpells x K
	popAtRisk.providerAdoptions = initialConditions.popAtRisk.providerAdoptions' - cumsum(providerDeltas, 1) + providerDeltas; % NumSpells x K
	popAtRisk.providerExits      = installedBases.laggedProviders;                               % NumSpells x K
	clear providerDeltas;
	
	% Make laggedConsumersOrig and laggedConsumersDest + popAtRisk for consumerAdoptions
	[Y_consumerAdoptionsOrig, Y_consumerAdoptionsDest] = aggregate_sparseY(Y_full.consumerAdoptions, NumSpells, K1, K2);
		% [NumSpells x K1] and [NumSpells x K2]
	installedBases.laggedConsumersOrig = initialConditions.laggedConsumersOrig' + cumsum(Y_consumerAdoptionsOrig, 1) - Y_consumerAdoptionsOrig; % NumSpells x K1
	installedBases.laggedConsumersDest = initialConditions.laggedConsumersDest' + cumsum(Y_consumerAdoptionsDest, 1) - Y_consumerAdoptionsDest; % NumSpells x K2		
	popAtRisk.consumerAdoptions = initialConditions.popAtRisk.consumerAdoptions' - cumsum(Y_consumerAdoptionsOrig, 1) + Y_consumerAdoptionsOrig; % NumSpells x K1
	clear Y_consumerAdoptionsOrig Y_consumerAdoptionsDest;

	
	% Reshape installedBases fields as vectors
	installedBases.laggedConsumersOrig = installedBases.laggedConsumersOrig(:); % (NumSpells*K1) x 1
	installedBases.laggedConsumersDest = installedBases.laggedConsumersDest(:); % (NumSpells*K2) x 1
	installedBases.laggedProviders      = installedBases.laggedProviders(:); % (NumSpells*K) x 1
	
	% Reshape Y_full elements as vectors
	Y_full.consumerAdoptions = Y_full.consumerAdoptions(:); % (NumSpells*K1*K2) x 1
	Y_full.providerAdoptions  = Y_full.providerAdoptions(:);  % (NumSpells*K) x 1
	Y_full.providerExits      = Y_full.providerExits(:);      % (NumSpells*K) x 1
	
	%%% Make object data based on data generated
	% Copy dims from data0
	data.consumerAdoptions.dims = data0.consumerAdoptions.dims;
	data.providerAdoptions.dims  = data0.providerAdoptions.dims;
	data.providerExits.dims      = data0.providerExits.dims;
	
	% Copy Xparts from data0
	data.consumerAdoptions.Xparts = data0.consumerAdoptions.Xparts;
	data.providerAdoptions.Xparts  = data0.providerAdoptions.Xparts;
	data.providerExits.Xparts      = data0.providerExits.Xparts;
	
	% Get Y's that were generated
	data.consumerAdoptions.Y = Y_full.consumerAdoptions;  % (NumSpells*K1*K2) x 1
	data.providerAdoptions.Y  = Y_full.providerAdoptions; % (NumSpells*K) x 1
	data.providerExits.Y      = Y_full.providerExits;     % (NumSpells*K) x 1
	
	%%%%% Process things for part on consumerAdoptions
	%% Expand data to ctsTime
	data.consumerAdoptions.dims.NumObs   = NumSpells*K1*K2;
	data.consumerAdoptions.dims.dims1{2} = NumSpells*K1;
	data.consumerAdoptions.dims.dims1{3} = NumSpells*K2;
	data.consumerAdoptions.dims.dims1{4} = NumSpells;
	data.consumerAdoptions.Xparts{4}.X     = data.consumerAdoptions.Xparts{4}.X(spell2period,:);
	data.consumerAdoptions.Xparts{4}.X_FEs = data.consumerAdoptions.Xparts{4}.X_FEs(spell2period,:);
	data.consumerAdoptions.spellDurations = spellDurations;
	data.consumerAdoptions.spell2period   = spell2period;
	
	%% Expand Xdynamic_staticInputs to ctsTime
	tmp2 = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X;
	tmp2 = reshape(tmp2, [T K1 size(tmp2,2)]);
	tmp2 = tmp2(spell2period,:,:);
	tmp2 = reshape(tmp2,[NumSpells*K1 size(tmp2,3)]);
	Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X = tmp2;
	clear tmp2;
	
	tmp2 = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X_FEs;
	tmp2 = reshape(tmp2, [T K1 size(tmp2,2)]);
	tmp2 = tmp2(spell2period,:,:);
	tmp2 = reshape(tmp2,[NumSpells*K1 size(tmp2,3)]);
	Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X_FEs = tmp2;
	clear tmp2;
	
	tmp3 = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X;
	tmp3 = reshape(tmp3, [T K2 size(tmp3,2)]);
	tmp3 = tmp3(spell2period,:,:);
	tmp3 = reshape(tmp3,[NumSpells*K2 size(tmp3,3)]);
	Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X = tmp3;
	clear tmp3;
	
	tmp3 = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X_FEs;
	tmp3 = reshape(tmp3, [T K2 size(tmp3,2)]);
	tmp3 = tmp3(spell2period,:,:);
	tmp3 = reshape(tmp3,[NumSpells*K2 size(tmp3,3)]);
	Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X_FEs = tmp3;
	clear tmp3;
	
	
	%%%%% Process things for part on providerAdoptions
	%% Expand data to ctsTime
	data.providerAdoptions.dims.NumObs   = NumSpells*K;
	data.providerAdoptions.dims.dims1{1} = NumSpells*K;
	data.providerAdoptions.dims.dims1{2} = NumSpells;
	data.providerAdoptions.Xparts{2}.X     = data.providerAdoptions.Xparts{2}.X(spell2period,:);
	data.providerAdoptions.Xparts{2}.X_FEs = data.providerAdoptions.Xparts{2}.X_FEs(spell2period,:);
	data.providerAdoptions.spellDurations = spellDurations;
	data.providerAdoptions.spell2period   = spell2period;
	
	%% Expand Xdynamic_staticInputs to ctsTime
	tmp1 = Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X;
	tmp1 = reshape(tmp1, [T K size(tmp1,2)]);
	tmp1 = tmp1(spell2period,:,:);
	tmp1 = reshape(tmp1,[NumSpells*K size(tmp1,3)]);
	Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X = tmp1;
	clear tmp1;
	
	
	%%%%% Process things for part on providerExits
	%% Expand data to ctsTime
	data.providerExits.dims.NumObs   = NumSpells*K;
	data.providerExits.dims.dims1{1} = NumSpells*K;
	data.providerExits.dims.dims1{2} = NumSpells;
	data.providerExits.Xparts{2}.X      = data.providerExits.Xparts{2}.X(spell2period,:);
	data.providerExits.Xparts{2}.X_FEs  = data.providerExits.Xparts{2}.X_FEs(spell2period,:);
	data.providerExits.spellDurations = spellDurations;
	data.providerExits.spell2period   = spell2period;
	
	%% Expand Xdynamic_staticInputs to ctsTime
	tmp1 = Xdynamic_staticInputs.providerExits.Xparts_static{1}.X;
	tmp1 = reshape(tmp1, [T K size(tmp1,2)]);
	tmp1 = tmp1(spell2period,:,:);
	tmp1 = reshape(tmp1,[NumSpells*K size(tmp1,3)]);
	Xdynamic_staticInputs.providerExits.Xparts_static{1}.X = tmp1;
	clear tmp1;
	
	
	%%%%% Reshape extraArgs (surface)
	tmp = Xdynamic_staticInputs.consumerAdoptions.extraArgs{1};
	tmp = reshape(tmp, [size(tmp,1)/K1 K1]);
	tmp = repmat(tmp(1,:), [NumSpells 1]); % NumSpells x K1
	tmp = reshape(tmp, [NumSpells*K1 1]); % (NumSpells*K1) x 1
	Xdynamic_staticInputs.consumerAdoptions.extraArgs{1} = tmp;
	Xdynamic_staticInputs.providerAdoptions.extraArgs{1}  = tmp;
	Xdynamic_staticInputs.providerExits.extraArgs{1}      = tmp;
	clear tmp;
	
	
	% Change the Xparts that are "dynamic", based on installed bases that were generated
	Xdynamics_consumerAdoptions = make_Xdynamic_consumerAdoptions(indivConfigs.consumerAdoptions, Xdynamic_staticInputs.consumerAdoptions, installedBases, true);
	data.consumerAdoptions.Xparts{2} = Xdynamics_consumerAdoptions{2};
	data.consumerAdoptions.Xparts{3} = Xdynamics_consumerAdoptions{3};
	clear Xdynamics_consumerAdoptions;
	
	Xdynamics_providerAdoptions = make_Xdynamic_providerAdoptions(indivConfigs.providerAdoptions, Xdynamic_staticInputs.providerAdoptions, installedBases, true);
	data.providerAdoptions.Xparts{1} = Xdynamics_providerAdoptions{1};
	clear Xdynamics_providerAdoptions;
	
	Xdynamics_providerExits = make_Xdynamic_providerExits(indivConfigs.providerExits, Xdynamic_staticInputs.providerExits, installedBases, true);
	data.providerExits.Xparts{1} = Xdynamics_providerExits{1};
	clear Xdynamics_providerExits;
	
	% Make logLambdaOffset
	laggedProviders = installedBases.laggedProviders; % (T*K) x 1
	data.providerAdoptions.logLambdaOffset = log(popAtRisk.providerAdoptions(:)); % (T*K) x 1
	data.providerExits.logLambdaOffset     = log(laggedProviders);             % (T*K) x 1
	data.providerExits.logLambdaOffset     = max(data.providerExits.logLambdaOffset, -2e4); % (T*K) x 1
	data.consumerAdoptions.logLambdaOffset_dim1     =  log(popAtRisk.consumerAdoptions); % T x K1
	data.consumerAdoptions.logLambdaOffset_dim2 = - 2e10 * reshape(laggedProviders <= 0, [data.consumerAdoptions.dims.dims1{4} 1 K2]); % T x 1 x K2
	
	% Make LL_offset for each submodel
	idxes = find(data.providerAdoptions.Y); 
	data.providerAdoptions.LL_offset = -sum(gammaln(1+data.providerAdoptions.Y(idxes)));
	clear idxes;
	
	idxes = find(data.providerExits.Y); 
	data.providerExits.LL_offset = -sum(gammaln(1+data.providerExits.Y(idxes)));
	clear idxes;

	idxes = find(data.consumerAdoptions.Y); 
	data.consumerAdoptions.LL_offset = -sum(gammaln(1+data.consumerAdoptions.Y(idxes)));
	clear idxes;
	
	%%% Save data to output file
	disp('Saving to output file...');
	rep_outputFolder = createFolderIfNotExist(sprintf('%s/replication_%03d', outputFolder, iter));
	outputMatFile = sprintf('%s/syntheticData.mat', rep_outputFolder);
	save(outputMatFile, 'data', 'truth', 'startingRandomState', '-v7.3');
	
	%%% Clean up
	clear Y_full installedBases laggedProviders popAtRisk;
	
	
	%%% Aggregate to day-level and save to second output file
	data.providerAdoptions  = aggregateCtsTimeData_providerAdoptions(data.providerAdoptions);
	data.providerExits      = aggregateCtsTimeData_providerExits(data.providerExits);
	data.consumerAdoptions  = aggregateCtsTimeData_consumerAdoptions(data.consumerAdoptions);
	
	rep_outputFolder2 = createFolderIfNotExist(sprintf('%s/replication_%03d', outputFolder2, iter));
	outputMatFile2 = sprintf('%s/syntheticData.mat', rep_outputFolder2);
	save(outputMatFile2, 'data', 'truth', 'startingRandomState', '-v7.3');
	
	% Clean up
	clear data;
end

end